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Dynamical mean field theory of optical third harmonic generation 
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We formulate the third harmonic generation (THG) within the dynamical mean field theory 
(DMFT) approximation of the Hubbard model. In the limit of large dimensions, where DMFT 
becomes exact, the vertex corrections to current vertices are identically zero, and hence the 
calculation of the THG spectrum reduces to a time-ordered convolution, foUowd by appropriate 
analytic continuuation. We present the typical THG spectrum of the Hubbard model obtained 
by this method. Within our DMFT calculation, we observe a nontrivial approximate scaling 
function describing the THG spectra in all Mott insulators, independent of the gap magnitude. 
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Nonlinear optical interactions of laser fields with mat- 
ter provide powerful spectroscopic tools for the under- 
standing of microscopic processes. The ability to con- 
trol pulse durations (to a few femtoseconds), bandwidths 
(up to 1 Hz resolution), and peak intensities (up to 
lO^^W/cm^) provides novel probes of elementary dy- 
namic events of matter.^ 

Observation of large third order nonlinear susceptibil- 
ity in a quasi one dimensional Mott insulator Sr2Cu03 



pulsion fall approximately on the same curve, if we scale 
the frequencies with the gap magnitude. This behavior 
is similar to the one observed in band insulators,'^ where 
a single particle picture describes the nonlinear optical 
processes. 

We start with the Hubbard model at half filling 
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(x(3) values in the range 10 » to 10 ^ e.s.u.)^'^ poses ^^^^^ ^t^ ^^g^tes an electron at site i with spin s 



the problem of nonlinear optical response in correlated 
insulators. 

In systems with large on-site Coulomb interaction, the 
ID system has the largest optical nonlinearity because of 
the decoupling of spin and charge degrees of freedom.^' ^ 
On the other hand, mean field analysis shows that among 
SDW-ordered systems, the largest third order optical re- 
sponse appears in 2D.^ A natural generalization to higher 
dimensions is through the dynamical mean field theory 
(DMFT). DMFT includes the effect of local quantum 
fluctuations and becomes exact in the limit of infinite 
dimensions. 

In this limit the self-energy becomes a purely lo- 
cal quantity determined from a self-consistent Ander- 
son model. Then the matrix elements determining spec- 
tral weights are encoded in the local self-energy. In this 
limit the matrix element summations reduce to density of 
states (DOS) integrations. Therefore the combined DOS 
and self-energy effects gives rise to nonlinear optical re- 
sponse of the system. 

We formulate a nonlinear response theory, with exam- 
ple of THG within DMFT approximation and prescribe 
a numerically feasible method to avoid expensive com- 
putations. To our knowledge this is the first application 
of DMFT to nonlinear optics. 

The general THG line shape within DMFT framework 
consists in a strong peak at three photon resonance, fol- 
lowed by a shoulder at two photon resonance, and a very 
weak feature at one photon resonance. The three-photon 
contributions obtained for various on-site Coulomb re- 
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The dimension of the lattice is d. Here the 1/Vd scaling 
of the hopping term ensures that average kineitc energy 
per particle in the limit of d — > oo remains finite.* 

If we now imagine that we integrate out all degrees 
of freedom on various lattice sites, except for the one at 
the origin, we will be left with an effective action for this 
"impurity" site: 

5eff - - / / drdr' J2 cUr)G-'{r - t')co.{t') 

^° ^° ^ (2) 

+ U \ not(r)nox(T). 
Jo 

Here the impurity propagator C/o(t— r') describes tempo- 
ral quantum fluctuations between the four possible states 
of a single site at the origin, which must be determined 
self-cosistently. 

The simplest way to solve such effective impurity prob- 
lem is the second order perturbation, which gives 



s(r) = ^g,(r)g,(T)g,{-T) 



(3) 



This gives the lattice Green's function as, 

G'(fc,ia;„) = l/(iw„ - £j: - i;(Ja;„)). (4) 

The projection of this function on site 'o' is given by 

D{z)de 



G{iLu„ 



(5) 



where D{e) is the lattice density of states. Finally the 
self-consistency between lattice (G) and impurity (Qo) is 
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via the Dyson equation 



(6) 



which is used to update Qo if the consistency has not 
been achieved yet.® 

Solving the set of equations (3), (5) and (6) for var- 
ious values of Hubbard U captures the physics of Mott 
metal-to-insulator transtion.^ The essential many-body 
quantity provided by solving the local impurity problem 
is the self-energy which encodes the matrix element ef- 
fects, as will be shown in the following. We solve^*' the 
above set of equations at zero temperature for a Bethe 
lattice DOS of type D{e) = -Vl — e^, which corresponds 

to renormalized hopping i — tyd = 1/2. 

The third order nonlinear optical response per 
unit volume is related to four-current correlation 

Xri{^;^ii^2,L^3) by" 



X '(uj;UJl,U}2,U}3) 
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^ 1 xf^{'^;^i,'^2,ujs) 



V 



tLJc,LJlUJ2UJ[i 



, (7) 



where lu = —uj„ = — (t^i + uj2 + W3), and the 4— current 
correlation function is given by 



X„- (w;(i;i,W2,W3) = / dxdxidx2dx^e 



iujt+u:it\+uj'2t2+oJzt2. 



X {Tc]{x)]{xi)j{x2)3{x3.)). 

where j{x) is current operator at space-time x = {r,t). 
Here Tc denotes the time-ordering along the Keldysh 
path.^^ Although the general formulation can be writ- 
ten down in terms of Keldysh Green's functions, but for 
parametric processes,^ i.e. processes in which, final and 
initial states are identical, in practice one can avoid use 
of Keldysh Green's functions. In such a case one can use 
ordinary Green's function, to calculate the time ordered 
diagrams, followed by appropriate analytic continuation 
to ensure the correct v + ir] behavior of the fully retarded 
optical responses.'^ 

The case of third harmonic generation corresponds to 
LOi ^ UJ2 = uj'i — V, so that we have 



X [v] = x ( 



-iv; v,v,v) = -— [ — 
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The Feynman diagram corresponding to the THG pro- 
cess^"^ is shown in Fig. 1. In the limit of infinite dimen- 
sions vertex corrections to odd parity operators indenti- 
cally vanish by Ward indentity. To see this, let us write 
down the Ward identity as^** 

-ik^Vip + k,k)^ G-^{p + k)- G-\k) (9) 

where summation over /i = 0, 1 . . . , d is understood and 
p, k are {d + l)-vectors. Using the Dyson equation the 
right hand side becomes T,(k+p) — S(fc). In d — > 00 limit 
self-energy is purely local® (no k dependence) and hence 
it vanishes. Now, since the current (ocvelocity) vertex 
has odd parity under k -^ —k, the vertex correction F^ 
identically vanishes. 

Therefore the four-current correlation {jjjj)'^^ {v) 
in Fig. 1 can be obtained by simple convolution. The 




Fig. 1. Feynman diagram corresponding to the third harmonic 
generation. This diagram represents the time-ordered four- 
current correlation. To obtain the fully retarded four-current cor- 
relation we ensure the correct v + ir) analytic behavior. 



Green's functions running around the loop are self- 
consistent lattice Green's function obtained from the im- 
purity solver by iterated perturbation theory.® Now let us 
further simplify {jjjj)'^^ {v) in d ^ cx) limit. Equation 
(8) can be written as 

1 (,7i.?i.7i.7i) 8F 



THG 
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— ^sin4(fci)x 



^fea-(*^")G'j^(itj„ + iv)G^^{iu!n + 2iv)G^^{iujn + iiv) 

where we have used the fact that current vertex in di- 
rection no. 1 is 2isin(/ci). To proceed further, we need 

to define po{£) = 1^ X^J si'^^('''i)'^(^ ~ ^i)- To take the 
limit d ^ 00, we Fourier transform po{£) as: 

Po(s) = J pois)e^''de = -i^^sin4(fci)e-^£ 
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where J„ is the Bessel function of order n. In the last line 
we have used the fact that in d — > cx) the hopping matrix 
element scales like t =■ tl\fd and hence st <^ 1, so that 
using J„(x) K, 2!— r, we can ignore J 2 and J4 compared 
to Jq. Repeating the above algebra without sin k\ shows 
that, [Jo(2si)] is indeed the Fourier transform of Dis)- 
Therefore 



Po(e) = o^(e)' 



which allows us to write 

?4 



X 



THG 



{^y) 



^SZ/'^^W 



(10) 



(11) 



G(e, iujn)G{£, iujn + iv)G{e, Jcj„ + 2iiy)G{e, Jcj„ + 3iiy) 

We see that in rf ^ 00 the k summation becomes simply 
a DOS integration. In the following e stands for e^, and 
the explicit e subscript emphasises the k label. From the 
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derivation it can be seen that other four current correla- 
tions like O1J1J2J2) in d ^ 00 limit differ from (jijijiji) 
by a numerical factor. Therefore to that extent the limit 
of d — » 00 is blind to various directions in space. Hence 
DMFT can not distinguish optical spectroscopics with 
polarized light from those with unpolarized light. 

To elucidate the matrix element effects in DMFT 
method, after Lehman representing the Greens' functions 
in terms of A{k, E) — > A{e^ E), and using standard con- 
tour integration techniques to perform l//3^j^ summa- 
tion we obtain 



+4 r 

X™°(iJ/) = -ir^ I deD{e)dEi ...dE^ 
A{e,Ei)A{e,E2)A{e,E^)A{e,Ei) x F, 



(12) 



with 



f{E^) 



El- E2+ iv El - 


- E3 + 2w El- Ei + 3iiy 
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1 1 



(13) 



E4-E1- iiv E4- E2- 2w E4- E3- iv 

where / is the Fermi function. This expression closely 
resembles familiar expressions in nonlinear optics litera- 
ture (see e.g. Sec. 3.2 of Ref. 2). Therefore in this for- 
mulation, the matrix element effects appear via spectral 
function A{e,E), which itself is fully determined by the 
self-energy. In principle after replacing iv with v + iQ^ in 
this expression, we can use the spectral weights obtained 
from DMFT solver to calculate the nonlinear response. 
However, numerical calculation of the above five dimen- 
sional integrals is not computationally feasible. 

Another alternative method would be to calculate 



X 



THG 



at Matsubara frequencies according to (11), fol- 
lowed by analytic continuation iv — > v + iQ'^ . But in this 
process, we face spurious features characteristic of ana- 
lytic continuation of numbers, which makes it difficult to 
assess the nonlinear dynamical structures. 

Since we are interested in high energy features in the 
scale of Mott gap, which is much larger than the thermal 
energies at room temperature, in order to avoid the above 
mentioned difficulties, we work at zero temperature. At 
T = 0, (11) will be replaced by 

,THG/ ^_ t^ f dudePje) 1 1 1 



X 



Gttv'^ 



where ^j — lu + jv — Sj^(cj + jv) + i|S/(cj + jv)\ for 
i = 0,1,2,3. 

In the above formula, (i) the integration over e cor- 
responds to summation over intermediate states in con- 
ventional expressions^ which are usually used for systems 
with discrete energy levels, and (ii) the matrix element 
effects are encoded in S(a;), real and imaginary part of 
which have been denoted by E/j and S/, respectively. It 
is very crucial to note that we have used absolute value 
of the imaginary part of the self-energy. This is indeed 
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Fig. 2. (Color online) THG and DOS for U = 4.5. We use the 
Bethe lattice in solving DMFT equations. In the upper pannel, 
dashed line is the imaginary part of (jjjj) , narrow solid line repre- 
sents its real part, and the bold solid line represents the absolute 
value. Lower pannel shows the DOS. The onset of (3-photon) 
absorption at u ^ 0.85 corresponds to the gap value, Eg sd 2.5, 
while peak-to-peak energy difference is on the scale of U = 4.5 



a necessary step to transform from time-ordered four- 
current, to fully retarded one.'^ 

Decomposing the integrand in (14) to partial fractions 
in terms of e, the four-current correlation can be written 
as /o + /i + /2 + /s, where 

/,(.) = £ /dc^nT^^ fdej^ (15) 



6tt 



ST^j 



(0-6)7 ie~Q 



The expression (15) reveals the resonance structure 
transparently: It arises from uj + jv — T,r{lli + jv) — e. 
When the frequency v of the photons is such that j- 
photon frequency matches the energy difference e — u + 
E/j(a; + jv), we will have a j-photon resonance. Here, /g 
corresponds to a background contribution. 

Note that within this formulation, we do not have any 
control over the 'broadening' parameter rj, instead the 
broadening 77 — |E/(cj)| is determined by the solution of 
the impurity problem in a self-consistent fashion. 

Now let us present our results. For semicircular Bethe 
lattice DOS of width 2i = 1, the critical value is given by 
Uc « 3.3, above which system is in the Mott insulating 
state. Fig. 2 shows the result for U = 4.5. Lower pannel 
shows the self-consistent DOS, with a Mott-Hubbard gap 
~ 2.5. The peak-to-peak separation between the upper 
and lower Hubbard bands is ^ U — 4.5. Upper pan- 
nel shows real(dashed), imaginary (dotted) and absolute 
value(solid line) of the four-current correlation {jjj])'^^'^ ■ 

The onset of absorption starts dX v k, 0.85 which 
is 1/3 of the gap magnitude. This can be clearly seen 
in the imaginary part of the THG four-current correla- 
tion in Fig. 2. This onset clearly corresponds to three- 
photon absorption. The three-photon resonance peaks 
around v k, \.h (denoted by A) which is 1/3 of the peak- 
to-peak separation of the Hubbard bands. Next weaker 
feature (denoted by B), which is a shoulder similar to 
finite-dimensional results,^ corresponds to 1/2 of peak- 
to-peak separation (^ 4.5) of the Hubbard bands. The 
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U=4.5 




200 
150 


1 'I 

'l I \> 

1 


A 

'A 


---l"^<iijCial 
lm<jjij>*9 


100 


A 

1 \ 


- 


50 




1 \ 

V \ 





1 








~^ 


-50 




\ ^ y^ 



In our previous work,'^ we found that within the SDW 
mean field the THG in e.g. ID is given by S5x"'"^*^(i^) = 



Fig. 3. Solid line shows the imaginary part of the total THG for 
U = 4.5. The dashed line shows the imaginary part of f3{v). A 
and B in this figure correspond to those in Fig. 2. 



one-photon process is the weakest feature around v ^ 4.5 
which can hardly be distinguished in THG spectrum. 
Despite that DMFT is designed to work better in larger 
dimensions, the spectrum in figure 2 qualitatively resem- 
bles the experimental result on Sr2Cu03, which is a one- 
dimensional Mott insulator.^ 

We claimed that in peak A of Fig. 2 the dominant con- 
tribution comes from three-photon processes. To demon- 
strate this, in Fig. 3 we plot with dashed line, the imag- 
inary part of f^{v). The solid line shows the total four- 
current correlation. As can be seen, the 3— photon reso- 
nance dominantly contributes to the peak A, although 
the peak position is slightly shifted to lower energies. 
Further, it is clearly seen that fz{v) deos not contribute 
much to the two photon resonance at B. /3(i^) also blows 
up at small frequencies which by (15) will be finally com- 
pensated by other patial spectra, /O, /I, /2, to give the 
total THG spectrum. 




f/(f)-8/M + i/(f), in which 
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Fig. 4. Imaginary part of fz{v) vs. u/ Eg, where Eg is the Mott- 
Hubbard gap in the Mott insulating phase. 



where the gap Eg = 2A, w'^ ^ {I + A2)/A2, and the 
frequency v = ^ is scaled by the gap parameter. In the 
above equation f{2>v/2) corresponds to fz{v) defined in 
(15). We see that in the mean field approximation, fz{j^) 
has a scaling part which is universal, and independent of 
the gap magnitude. 

Motivated by this observation, in Fig. 4 we plot Im/3 
for Mott insulators with various values of C/ as a func- 
tion of -^. We also apply an overal scaling to the curves. 
Such a scaling behavior, though approximate, points to a 
universal features in the nonlinear optical spectra of the 
Mott insulators, which are independent of the gap mag- 
nitude. It seems that the mean field scaling behavior sur- 
vives the quantum fluctuations implemented via DMFT. 
It would be interesting to further explore this observa- 
tion using alternative methods of dealing with correlated 
insulators. 

In summary, on the technical side, we have formulated 
the nonlinear optical response theory within the DMFT 
theory. In equation (14) we present a feasible way that 
avoids numerical difficulties. On the physics side, assum- 
ing that DMFT is a good approximation for d > 1, we 
observe that nonlinear optical spectra of higher dimen- 
sional Mott insulators share common features with those 
observed in rf = 1 dimensional Mott insulator, Sr2Cu03.^ 
Also within DMFT we observe an approximate scaling 
behavior in the THG spectra of Mott insulators. 
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